This script carries out the beta diversity analyses (community composition and structure) for the 2016 Cork great tit microbiome data.
Overall plan" * model 1: all variables of interest, nest as blocking factor * model 2: all variables of interest but age-habitat interaction, nest blocking * model 3: adult samples, all variables of interest, nest blocking * then check are results different when individuals w/ 2 samples are dropped
Gabrielle says not to remove individuals sampled twice
library(devtools)
Loading required package: usethis
Attaching package: ‘devtools’
The following object is masked from ‘package:permute’:
check
install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
Skipping install of 'pairwiseAdonis' from a github remote, the SHA1 (ece560d2) has not changed since last install.
Use `force = TRUE` to force installation
phyloseq-class experiment-level object
otu_table() OTU Table: [ 54343 taxa and 246 samples ]
sample_data() Sample Data: [ 246 samples by 50 sample variables ]
tax_table() Taxonomy Table: [ 54343 taxa by 7 taxonomic ranks ]
refseq() DNAStringSet: [ 54343 reference sequences ]
Remove individuals measured twice * set.seed so same individuals dropped * no adults measured twice
#n_occur <- data.frame(table(metadata.NoDups$bird.ID))
#n_occur[n_occur$Freq > 1,]
set.seed(1189)
metadata.NoDups <- metadata %>%
group_by(bird.ID) %>%
sample_n(1)
dropped.samples <- setdiff(metadata$BIOM.ID, metadata.NoDups$BIOM.ID)
# metadata.even <- metadata %>%
# group_by(bird.ID) %>%
# sample_n(min(table(metadata$bird.ID)))
table(metadata$ageBinned)
1week 2week adult
81 114 51
table(metadata.NoDups$ageBinned)
1week 2week adult
55 91 51
Do i need to scale numeric variables?
numeric.predictors <- c("QubitDNA","Tarsus","Weight", "wing", "broodSizeWhenSampled", "broodSizeMax", "totalFledge", "clutchSize", "layDateFirst", "numberDeadPreRinged", "numberDeadPostRinged","scaled.mass","scaled.mass.wing.adult", "scaled.mass.tarsus.adult","scaled.mass.chick", "scaled.mass.tarsus", "scaled.mass.wing", "DistanceToEdge")
# centre and scale numeric variables ie. subtract mean and divide by st. deviation
metadata.NoDups.scaled <- metadata.NoDups
Error: object 'metadata.NoDups' not found
Compositional method, use clr transform
Keeping repeated measures here but using blocking factor to control for repeated samples.
# Filter rare taxa
phylo.knowles <- filter_taxa(phylo.spring, function(x) sum(x > 1) > (0.05*length(x)), TRUE)
# phylo.knowles.clr <- phylo.knowles
# phylo.knowles.clr@otu_table <- otu_table(clr(phylo.knowles@otu_table), taxa_are_rows = F)
# aitchison.dist.dups <- phyloseq::distance(phylo.knowles.clr, method = "euclidean")
pk.otu.clr.dups <- clr(phylo.knowles@otu_table)
aitchison.dist.dups <- vegdist(pk.otu.clr.dups, method = "euclid")
# phylo.TSS <- transform_sample_counts(phylo.knowles, function(x) x/sum(x)) # normalise read counts w/ Total-Sum Scaling
#
# BCdist <- phyloseq::distance(phylo.TSS, method="bray")
All variables except for habitat have heterogenous dispersions
## H0= No difference in dispersion between groups
# calc dispersion, using distance measure
variables <- c("ageBinned", "habitat", "layDateFirst", "broodSizeWhenSampled", "DistanceToEdge", "SequencePlate")
for(i in variables){ # works
dispersion <- betadisper(aitchison.dist.dups, metadata.scaled[,i]) #, bias.adjust = T
print(i) # print variable being tested
print(permutest(dispersion, pairwise=FALSE, permutations=1000))
cat("\n") # print line break, makes it easier to read
}
#hist(metadata$layDateFirst)
#ist(sqrt(metadata$DistanceToEdge))
dispersion1 <- betadisper(aitchison.dist.dups, metadata$ageBinned) #, bias.adjust = T
print(permutest(dispersion1, pairwise=FALSE, permutations=1000))
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 107.0 53.5 3.7152 1000 0.02498 *
Residuals 243 3499.2 14.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
dispersion2 <- betadisper(aitchison.dist.dups, metadata$habitat) #, bias.adjust = T
print(permutest(dispersion2, pairwise=FALSE, permutations=1000))
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 5.7 5.7331 0.3846 1000 0.5504
Residuals 244 3637.6 14.9082
dispersion3 <- betadisper(aitchison.dist.dups, (metadata$layDateFirst)) #, bias.adjust = T
print(permutest(dispersion3, pairwise=FALSE, permutations=1000))
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 27 1047.2 38.785 2.6146 1000 0.000999 ***
Residuals 218 3233.8 14.834
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
dispersion4 <- betadisper(aitchison.dist.dups, metadata$broodSizeWhenSampled) #, bias.adjust = T
print(permutest(dispersion4, pairwise=FALSE, permutations=1000))
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 6 220.0 36.674 2.505 1000 0.02498 *
Residuals 239 3499.1 14.641
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
dispersion5 <- betadisper(aitchison.dist.dups, (metadata$DistanceToEdge)) #, bias.adjust = T
print(permutest(dispersion5, pairwise=FALSE, permutations=1000))
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 60 2515.5 41.924 4.0906 1000 0.000999 ***
Residuals 185 1896.0 10.249
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
dispersion6 <- betadisper(aitchison.dist.dups, metadata$SequencePlate) #, bias.adjust = T
print(permutest(dispersion6, pairwise=FALSE, permutations=1000))
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 4 541.68 135.419 10.509 1000 0.000999 ***
Residuals 241 3105.49 12.886
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(dispersion1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = distances ~ group, data = df)
$group
diff lwr upr p adj
2week-1week -1.4854366 -2.785847 -0.1850262 0.0205827
adult-1week -1.1177060 -2.717328 0.4819163 0.2277401
adult-2week 0.3677306 -1.139788 1.8752487 0.8334642
TukeyHSD(dispersion2)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = distances ~ group, data = df)
$group
diff lwr upr p adj
deciduous-conifer 0.3765808 -0.8195671 1.572729 0.5357527
#TukeyHSD(dispersion3)
#TukeyHSD(dispersion4)
#TukeyHSD(dispersion5)
Plot dispersions
perms.dups <- with(metadata, how(nperm = 1000, blocks = nest))
all.adonis.dups.fixed <- adonis2(aitchison.dist.dups ~ ageBinned + habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, by="margin", method="euclidian", data = metadata.scaled, permutations = perms.dups)
pairwiseAdonis::pairwise.adonis2(aitchison.dist.dups ~ ageBinned + habitat, data = metadata.scaled)
$parent_call
[1] "aitchison.dist.dups ~ ageBinned + habitat , strata = Null , permutations 999"
$`2week_vs_adult`
Df SumOfSqs R2 F Pr(>F)
ageBinned 1 647 0.01699 2.8293 0.001 ***
habitat 1 395 0.01036 1.7263 0.001 ***
Residual 162 37034 0.97265
Total 164 38075 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
$`2week_vs_1week`
Df SumOfSqs R2 F Pr(>F)
ageBinned 1 635 0.01317 2.5796 0.001 ***
habitat 1 335 0.00695 1.3625 0.027 *
Residual 192 47228 0.97988
Total 194 48198 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
$adult_vs_1week
Df SumOfSqs R2 F Pr(>F)
ageBinned 1 563 0.01645 2.1756 0.001 ***
habitat 1 270 0.00788 1.0425 0.365
Residual 129 33369 0.97566
Total 131 34201 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
attr(,"class")
[1] "pwadstrata" "list"
#?pairwise.adonis2()
Try same as above but strata by sequence plate and then inlcude individual ID
# perms.dups <- with(metadata, how(nperm = 1000, blocks = SequencePlate))
#
# all.adonis.dups.fixed <- adonis2(aitchison.dist.dups ~ ageBinned + habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + bird.ID, by="margin", method="euclidian", data = metadata.scaled, permutations = perms.dups)
#
# all.adonis.dups.fixed
Including bird.ID as a fixed effect and blocking by nest suggests age is in fact a significant factor though accounts for only 0.7% of variation, while bird id accounts for 75%, but is non-significant. Sequence plate is significant, accounting for 2.5%, brood size is not significant while habitat, lay date and distance to edge could not be estimated.
# perms.dups <- with(metadata, how(nperm = 1000, blocks = nest))
#
# all.adonis.dups.fixed <- adonis2(aitchison.dist.dups ~ ageBinned + habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate + bird.ID, by="margin", method="euclidian", data = metadata.scaled, permutations = perms.dups)
#
# all.adonis.dups.fixed
# perms.dups <- with(metadata, how(nperm = 1000, blocks = bird.ID))
#
# all.adonis.dups.fixed <- adonis2(aitchison.dist.dups ~ ageBinned + habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, by="margin", method="euclidian", data = metadata.scaled, permutations = perms.dups)
#
# all.adonis.dups.fixed
perms.dups <- with(metadata, how(nperm = 1000, blocks = nest))
all.adonis.dups.int <- adonis2(aitchison.dist.dups ~ ageBinned*habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, by="margin", data = metadata.scaled, permutations = perms.dups) #
all.adonis.dups.int
Permutation test for adonis under reduced model
Marginal effects of terms
Blocks: nest
Permutation: free
Number of permutations: 1000
adonis2(formula = aitchison.dist.dups ~ ageBinned * habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, data = metadata.scaled, permutations = perms.dups, by = "margin")
Df SumOfSqs R2 F Pr(>F)
layDateFirst 1 420 0.00693 1.8300 0.089910 .
broodSizeWhenSampled 1 305 0.00504 1.3312 0.016983 *
DistanceToEdge 1 502 0.00830 2.1901 0.000999 ***
SequencePlate 4 3634 0.06002 3.9609 0.000999 ***
ageBinned:habitat 2 595 0.00982 1.2960 0.010989 *
Residual 233 53444 0.88259
Total 245 60554 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Interaction calc by term rather than margin to get estimates for main effects as well as interaction term
perms.dups <- with(metadata, how(nperm = 1000, blocks = nest))
adonis2(aitchison.dist.dups ~ ageBinned*habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, by="term", data = metadata.scaled, permutations = perms.dups) #
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks: nest
Permutation: free
Number of permutations: 1000
adonis2(formula = aitchison.dist.dups ~ ageBinned * habitat + layDateFirst + broodSizeWhenSampled + DistanceToEdge + SequencePlate, data = metadata.scaled, permutations = perms.dups, by = "term")
Df SumOfSqs R2 F Pr(>F)
ageBinned 2 1239 0.02046 2.7003 0.000999 ***
habitat 1 347 0.00573 1.5124 0.000999 ***
layDateFirst 1 490 0.00810 2.1377 0.000999 ***
broodSizeWhenSampled 1 311 0.00514 1.3578 0.078921 .
DistanceToEdge 1 475 0.00784 2.0689 0.000999 ***
SequencePlate 4 3653 0.06033 3.9816 0.000999 ***
ageBinned:habitat 2 595 0.00982 1.2960 0.009990 **
Residual 233 53444 0.88259
Total 245 60554 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
table(metadata$ageBinned, metadata$habitat)
conifer deciduous
1week 16 65
2week 27 87
adult 8 43
For loop sample equal numbers of from each habitat and age combination
aitch.pca <- prcomp(aitchison.dist.dups)
#plot(aitch.pca)
biplot(aitch.pca)
autoplot(aitch.pca, data = metadata)#, colour = "ageBinned", loadings = T, loadings.label = T, scale = 0)
Code from: https://huboqiang.cn/2016/03/03/RscatterPlotPCA
df_out <- as.data.frame(aitch.pca$x)
df_out$ageBinned <- metadata$ageBinned
df_out$habitat <- metadata$habitat
df_out$layDateFirst <- metadata$layDateFirst
df_out$broodSizeWhenSampled <- metadata$broodSizeWhenSampled
head(df_out)
Calculate the PC loadings for labels * took these from different sites and they disagree alot
# Sum the total variance
d.mvar <- sum(aitch.pca$sdev^2)
# Calculate the PC1 and PC2 variance
PC1.label <- paste("PC1: ", round(sum(aitch.pca$sdev[1]^2)/d.mvar, 3)*100,"%")
PC2.label <- paste("PC2: ", round(sum(aitch.pca$sdev[2]^2)/d.mvar, 3)*100,"%")
#percentage <- round(aitch.pca$sdev / sum(aitch.pca$sdev) * 100, 2)
#percentage <- paste( colnames(df_out), "(", paste( as.character(percentage), "%", ")", sep="") )
#screeplot(aitch.pca)
Plot * Use same name and label for fill and shape in order to overlap legends * ellipse default is 95% confidence level for t distribution * can draw euclid ellipse where level = radius, but how to choose radius? * alpha specifies the transparency * https://stats.stackexchange.com/questions/217374/real-meaning-of-confidence-ellipse
p <- ggplot(df_out,aes(x=PC1,y=PC2, color=ageBinned ))
p <- p + geom_point(aes(shape = habitat)) + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1,
aes(group = interaction(ageBinned, habitat), fill = habitat, colour = ageBinned))
p <- p + scale_fill_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_shape_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_color_discrete(name = "Age", labels = c("Day 8", "Day 15", "Adult"))
#p <- p + scale_shape_discrete(element_blank())#, labels = NULL)
#p + facet_grid(habitat~.)
#theme<-theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"))
p+theme_classic()
Try euclidian ellipses, specify radius with level
p <- ggplot(df_out,aes(x=PC1,y=PC2, color=ageBinned ))
p <- p + geom_point(aes(shape = habitat)) + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="euclid", alpha=0.1,level = 10,
aes(group = interaction(ageBinned, habitat), fill = habitat, colour = ageBinned))
p <- p + scale_fill_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_shape_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_color_discrete(name = "Age", labels = c("Day 8", "Day 15", "Adult"))
#p <- p + scale_shape_discrete(element_blank())#, labels = NULL)
#p + facet_grid(habitat~.)
#theme<-theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"))
p+theme_classic()
facet.labels <- c("Conifer", "Deciduous")
names(facet.labels) <- c("conifer", "deciduous")
p + theme_classic() + facet_grid(habitat~., labeller = labeller(habitat = facet.labels))
# plot_ordination(phylo.TSS, ordBC, color = "ageBinned") +
# geom_point(size=3) +
# ggtitle("Ord: Bray-Curtis") + stat_ellipse(geom = "polygon", type="t", alpha=0.1, aes(fill=ageBinned))+
# theme_bw() + theme(panel.grid.major = element_blank(),
# panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
# panel.border = element_rect(linetype = "solid", colour = "black", size=.8)) +
# theme(text=element_text(size=14, family="serif"), axis.ticks = element_line(colour = "black", size = .7))+
# geom_point()
Plot features that contribute to classification
df_out_r <- as.data.frame(aitch.pca$rotation)
df_out_r$feature <- row.names(df_out_r)
df_out_r$ageBinned <- metadata$ageBinned
df_out_r$habitat <- metadata$habitat
df_out_r$layDateFirst <- metadata$layDateFirst
df_out_r$broodSizeWhenSampled <- metadata$broodSizeWhenSampled
df_out_r$DistanceToEdge <- metadata$DistanceToEdge
df_out_r
p<-ggplot(df_out_r,aes(x=PC1,y=PC2,label=ageBinned,color=habitat ))
#p<-p+geom_point()+theme + geom_text(size=3) + theme(legend.position = "none")
p
Plot one plot for each habitat ((ellipses each age) * change aplha to 0.05?
p <- ggplot(df_out,aes(x=PC1,y=PC2, color=ageBinned ))
p <- p + geom_point() + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1,
aes(group = interaction(ageBinned), fill = ageBinned, colour = ageBinned))
p <- p + scale_color_discrete(name = "Age", labels = c("Day 8", "Day 15", "Adult"))
p <- p + scale_fill_discrete(element_blank(), labels = NULL)
#p + facet_grid(habitat~.)
p+theme_classic()
facet.labels <- c("Conifer", "Deciduous")
names(facet.labels) <- c("conifer", "deciduous")
p + theme_classic() + facet_grid(habitat~., labeller = labeller(habitat = facet.labels))
one plot for each age (ellipses each habitat)
p <- ggplot(df_out,aes(x=PC1,y=PC2, color=habitat ))
p <- p + geom_point() + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1,
aes(group = (habitat), fill = habitat, colour = habitat))
p <- p + scale_color_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_fill_discrete(element_blank(), labels = NULL)
#p + facet_grid(habitat~.)
p+theme_classic()
facet.labels <- c("Day 8", "Day 15", "Adult")
names(facet.labels) <- c("1week", "2week", "adult")
p + theme_classic() + facet_grid(ageBinned~., labeller = labeller(ageBinned = facet.labels))
one plot for each age (ellipses each habitat)
p <- ggplot(df_out,aes(x=PC1,y=PC2, color=ageBinned ))
p <- p + geom_point() + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1,
aes(group = (ageBinned), fill = ageBinned, colour = ageBinned))
p <- p + scale_color_discrete(name = "Age", labels = c("D8", "D15", "Adult"))
p <- p + scale_fill_discrete(element_blank(), labels = NULL)
#p + facet_grid(habitat~.)
p+theme_classic()
p <- ggplot(df_out,aes(x=PC1,y=PC2, color=habitat ))
p <- p + geom_point() + xlab(PC1.label) + ylab(PC2.label) #, show.legend = F
p <- p + stat_ellipse(geom = "polygon", type="t", alpha=0.1,
aes(group = (habitat), fill = habitat, colour = habitat))
p <- p + scale_color_discrete(name = "Habitat", labels = c("Conifer", "Deciduous"))
p <- p + scale_fill_discrete(element_blank(), labels = NULL)
#p + facet_grid(habitat~.)
p+theme_classic()
p.habitat <- p
theme_bw(base_size = base_size, base_family = base_family) %+replace%
theme(
legend.position = "none",
panel.border = element_blank(),
axis.line = element_line(colour = "black"),
panel.grid.major=element_line(colour="grey", size=0.5, 3),
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
strip.background = element_blank(),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.5, face = "bold") # centre and bold title
)
Error in theme_grey(base_size = base_size, base_family = base_family, :
object 'base_size' not found
envfit1 <- envfit(aitch.pca,metadata$broodSizeWhenSampled, main="", col="red")
#plot(aitch.pca, type = "n")
plot(envfit1)
Error in strwidth(labels, ...) : plot.new has not been called yet
# from Jennys email + https://chrischizinski.github.io/rstats/ordisurf/
species.scores <- as.data.frame(scores(aitch.pca, "species")) # have to keep species here as the taxa
Error in names(x) <- paste("V", seq_along(x), sep = "") :
'names' attribute [1] must be the same length as the vector [0]
https://userweb.eng.gla.ac.uk/umer.ijaz/bioinformatics/ecological.html
#p<-p+theme(legend.key = element_blank(), #removes the box around each legend item
# legend.position = "bottom", #legend at the bottom
# legend.direction = "horizontal",
# legend.box = "horizontal",
# legend.box.just = "centre")
p
Error in rep(x[[i]], n) :
attempt to replicate an object of type 'closure'
Aitchison distance is supposedly robust to subsetting- can i therefore just subset the adult birds from previous distance calculation?
phylo.kn.adult <- filter_taxa(phylo.adults, function(x) sum(x > 1) > (0.05*length(x)), TRUE)
# phylo.TSS.adults <- transform_sample_counts(phylo.kn.adult, function(x) x/sum(x))
#
# BCdist.adults <- phyloseq::distance(phylo.TSS.adults, method="bray")
# JDdist.adults <- phyloseq::distance(phylo.TSS.adults, method="jaccard")
pk.otu.clr.adult <- clr(phylo.kn.adult@otu_table)
aitchison.dist.adult <- vegdist(pk.otu.clr.adult, method = "euclidean")
Homogenous dispersion: sex, habitat, distance is almost non-homo Non-homogenous: ageDays, …all others
## H0= No difference in dispersion between groups
# calc dispersion, using distance measure
variables <- c("ageDays", "Sex", "habitat", "layDateFirst", "broodSizeWhenSampled", "DistanceToEdge", "SequencePlate")
for(i in variables){ # works
dispersion <- betadisper(aitchison.dist.adult, adults.meta.scaled[,i])
print(i) # print variable being tested
print(permutest(dispersion, pairwise=FALSE, permutations=1000))
cat("\n") # print line break, makes it easier to read
}
Sex * habitat * broodSizeWhenSampled + Sex * habitat * layDateFirst + habitat * DistanceToEdge + (1|nest) + (1|SequencePlate)
perms.adult <- with(adults.meta, how(nperm = 1000, blocks = nest))
adult.adonis <- adonis2(aitchison.dist.adult ~ ageDays + Sex + habitat + DistanceToEdge + layDateFirst + broodSizeWhenSampled + SequencePlate, by = "margin", data = adults.meta.scaled, permutations = perms.adult)
adult.adonis
Permutation test for adonis under reduced model
Marginal effects of terms
Blocks: nest
Permutation: free
Number of permutations: 1000
adonis2(formula = aitchison.dist.adult ~ ageDays + Sex + habitat + DistanceToEdge + layDateFirst + broodSizeWhenSampled + SequencePlate, data = adults.meta.scaled, permutations = perms.adult, by = "margin")
Df SumOfSqs R2 F Pr(>F)
ageDays 1 271.0 0.01968 1.0475 0.629371
Sex 1 213.0 0.01546 0.8231 0.844156
habitat 1 290.3 0.02108 1.1221 0.449550
DistanceToEdge 1 327.0 0.02375 1.2640 0.315684
layDateFirst 1 323.6 0.02350 1.2507 0.145854
broodSizeWhenSampled 1 223.3 0.01621 0.8631 0.001998 **
SequencePlate 3 1462.1 0.10616 1.8836 0.003996 **
Residual 41 10608.1 0.77026
Total 50 13772.1 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#adonis2(JDdist.adults ~ Sex + habitat + DistanceToEdge + layDateFirst, by = "margin", data = adults.meta, permutations = perms.adult)
Table of all results for fixed and interaction PERMANOVA plus adult
#make df of lmer output
all.bdiv <- as_tibble(rbind(all.adonis.dups.fixed, all.adonis.dups.int, adult.adonis), rownames="Independent variables") %>%
dplyr::rename("P_estimate"="Pr(>F)") %>%
mutate_if(is.numeric, round, 3) %>%
mutate(P_estimate=ifelse(P_estimate==0,"<0.001",P_estimate)) %>%
mutate(P_estimate=ifelse(P_estimate<=0.05,str_c(P_estimate," *"),P_estimate)) %>%
mutate(P_estimate=ifelse(P_estimate>=0.05 & P_estimate<=0.06,str_c(P_estimate," ."),P_estimate))
# make df into kable
## this is in html, dosnt render in word doc # but can copy-paste
all.bdiv.results <- kable(all.bdiv, format = "html", table.attr = "style = \"color: black;\"") %>%
kableExtra::kable_styling(full_width = F) #%>%
# kableExtra::group_rows("(a) Top model, proteobacteria",1,nrow(proteobacteria.modelAvg.df)) %>%
#
#save_kable("alpha-kable__________.png") # this line saves as .png in Reports/
all.bdiv.results